Sea ice thickness

In this notebook, we’ll show more detailed analysis and visualization of ICESat2 sea ice thickness measurements and the PIOMAS (Pan-Arctic Ice Ocean Modeling and Assimilation System) modelled sea ice thickness data.

import xarray as xr # For working with gridded climate data 
import pandas as pd
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket

# Plotting dependencies
import cartopy.crs as ccrs
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore') 

1) Read in the data

We’ve included a function in the read_data_utils module that reads in the book dataset from the public google storage bucketas an xarray.Dataset object.

print(read_book_data.__doc__) # Print docstring
 Read in data for ICESat2 jupyter book. 
    If the file does not already exist on the user's local drive, it is downloaded from the books google storage bucket (https://console.cloud.google.com/storage/browser/is2-pso-seaice)
    The netcdf file is then read in as an xr.Dataset object 
    
    Args: 
        None
    Returns: 
        book_ds (xr.Dataset): data 
    
    
book_ds = read_book_data() # Read/download the data 
print(book_ds) # Show dataset
<xarray.Dataset>
Dimensions:                  (time: 30, y: 448, x: 304)
Coordinates:
  * time                     (time) datetime64[ns] 2018-11-01 ... 2021-04-01
    longitude                (y, x) float32 ...
    latitude                 (y, x) float32 ...
    region_mask              (y, x) uint8 ...
Dimensions without coordinates: y, x
Data variables: (12/18)
    freeboard                (time, y, x) float32 ...
    ice_density              (time, y, x) float32 ...
    ice_thickness            (time, y, x) float32 ...
    ice_thickness_smoothed   (time, y, x) float64 ...
    ice_thickness_unc        (time, y, x) float32 ...
    ice_type                 (time, y, x) float32 ...
    ...                       ...
    piomas_ice_thickness     (time, y, x) float32 ...
    t2m                      (time, y, x) float32 ...
    msdwlwrf                 (time, y, x) float32 ...
    drifts_uT                (time, y, x) float32 ...
    drifts_vT                (time, y, x) float32 ...
    drifts_magnitude         (time, y, x) float32 ...
Attributes:
    description:    Data used in ICESat-2 jupyter book
    note:           See individual data variables for references
    creation date:  2021-09-20

2) ICESat-2 interpolated sea ice thickness maps

We’ve used a linear interpolation method to generate a smoothed sea ice thickness variable and included it in the book’s dataset. Here, we’ll plot mean sea ice thickness for the winter seasons defined by Nov 2018 - Apr 2019, Sep 2019 - Apr 2020, and Sep 2020 - Apr 2021

def interactive_map(da, clabel="", cmap="viridis", colorbar=True, vmin=0, vmax=4, title="", frame_width=250): 
    pl = da.hvplot.quadmesh(y="latitude", x="longitude",
                            projection=ccrs.NorthPolarStereo(central_longitude=-45), 
                            features=["coastline"], # Add coastlines 
                            colorbar=colorbar, clim=(vmin,vmax), cmap=cmap, clabel=clabel, # Colorbar settings 
                            project=True, ylim=(60,90), 
                            dynamic=False, # Load the data before displaying
                            title=title, frame_width=frame_width) 
    return pl 
var = "ice_thickness_smoothed"

winter18_19 = pd.date_range(start="Nov 2018", end="Apr 2019", freq="MS")
winter19_20 = pd.date_range(start="Sep 2019", end="Apr 2020", freq="MS")
winter20_21 = pd.date_range(start="Sep 2020", end="Apr 2021", freq="MS")

mean18_19 = book_ds[var].sel(time=winter18_19).mean(dim="time",keep_attrs=True)
mean19_20 = book_ds[var].sel(time=winter19_20).mean(dim="time",keep_attrs=True)
mean20_21 = book_ds[var].sel(time=winter20_21).mean(dim="time",keep_attrs=True)

var_longname=book_ds[var].attrs["long_name"]
pl18_19 = interactive_map(mean18_19, vmin=0, vmax=4, colorbar=False, clabel=var_longname, frame_width=270, title="Mean Nov 2018 - Apr 2019") 

pl19_20 = interactive_map(mean19_20, vmin=0, vmax=4, colorbar=False, clabel=var_longname, frame_width=270, title="Mean Sep 2019 - Apr 2020") 

pl20_21 = interactive_map(mean20_21, vmin=0, vmax=4, clabel=var_longname, frame_width=270, title="Mean Sep 2020 - Apr 2021") 

data_pl = (pl18_19+pl19_20+pl20_21).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
data_pl

3) PIOMAS differences in winter sea ice thickness

Here we will compare the PIOMAS modelled sea ice thickness across winter seasons.

var = "piomas_ice_thickness"
var_longname=book_ds[var].attrs["long_name"]
book_ds[var] = book_ds[var].where(book_ds[var]>0)

mean18_19 = book_ds[var].sel(time=winter18_19).mean(dim="time",keep_attrs=True)
mean19_20 = book_ds[var].sel(time=winter19_20).mean(dim="time",keep_attrs=True)
mean20_21 = book_ds[var].sel(time=winter20_21).mean(dim="time",keep_attrs=True)
pl18_19 = interactive_map(mean18_19, vmin=0, vmax=4, colorbar=False, clabel=var_longname, title="Mean Nov 2018 - Apr 2019") 
pl19_20 = interactive_map(mean19_20, vmin=0, vmax=4, clabel=var_longname, title="Mean Sep 2019 - Apr 2020") 

diff = mean19_20 - mean18_19
diff_pl = interactive_map(diff, vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference") 

data_pl = (pl18_19+pl19_20+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl
pl19_20 = interactive_map(mean19_20, vmin=0, vmax=4, colorbar=False, clabel=var_longname, title="Mean Sep 2019 - Apr 2020") 
pl20_21 = interactive_map(mean20_21, vmin=0, vmax=4, clabel=var_longname, title="Mean Sep 2020 - Apr 2021") 

diff = mean20_21 - mean19_20
diff_pl = interactive_map(diff, vmin=-1.5, vmax=1.5, clabel="gridcell difference", cmap="coolwarm", title="Gridcell difference") 

data_pl = (pl19_20+pl20_21+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl

4) Mean monthly sea ice thickness comparison

Compare PIOMAS sea ice thickness and ICESat2 sea ice thickness

mean_monthly_pio = book_ds["piomas_ice_thickness"].mean(dim=["x","y"],keep_attrs=True)
mean_monthly_is2 = book_ds["ice_thickness_smoothed"].mean(dim=["x","y"],keep_attrs=True)
pio_lineplot = mean_monthly_pio.hvplot.line(grid=True, legend='top_left', ylabel="sea ice thickness") * mean_monthly_pio.hvplot.scatter(marker='o') # Overlay scatter plot to add markers
is2_lineplot = mean_monthly_is2.hvplot.line(grid=True, legend='top_left', ylabel="sea ice thickness") * mean_monthly_is2.hvplot.scatter(marker='o') # Overlay scatter plot to add markers
comb_plot = (pio_lineplot*is2_lineplot).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot